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We examine the correlation between compact radio quasars (redshifts in the range 
z = 0.3 — 2.2) and the arrival direction of ultrahigh energy cosmic rays forming clusters. 
Our Monte Carlo simulation reveals a statistically significant correlation on the AG ASA 
sample: the chance probability of this effect being less than 1%. The implications of 
this result on the origin of ultrahigh energy cosmic rays are discussed. 
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1 Introduction 



To date, some 20 giant air showers have been detected confirming the arrival of cosmic 
rays (CRs) with nominal energies at or above 10 20 ± 30% eV, with the record Fly's 
Eye event having ~ 3.2 x 10 20 eV Q. The mechanism(s) responsible for endowing 
particles with such enormous energies continues to present a major enigma to high 
energy physics 0. As shown in the pioneering works of Greisen-Zatsepin-Kuz'min 
(GZK) ||, the possible sources are constrained by the observed particle spectra due 
to the interaction with the universal radiation and magnetic fields on the way to the 
observer. In particular, any proton energy above 5 x 10 19 eV is degraded by resonant 
scattering via j+p — > A — ► p/n + n, such that less than 20% of protons survive with an 
energy above 3 x 10 20 (1 x 10 20 ) eV for a distance of 18 (60) Mpc. A typical nucleus of 
the cosmic radiation is subject to photodisintegration from blue-shifted relic photons, 
losing about 3-4 nucleons per traveled Mpc. The mean free path of gamma rays on the 
radio background decreases even more readily. Therefore, if the CR sources are all at 
cosmological distances, the observed spectrum must virtually end with the GZK cutoff 
at E « 8 x 10 19 eV. The spectral cutoff is less sharp for nearby sources (within 50 Mpc 
or so). The arrival directions of the trans-GZK events are distributed widely over the 
sky, without apparent counterparts (such as sources in the Galactic Plane or in the 
Local Supercluster) . Moreover, the data are consistent with an isotropic distribution 
of sources, in sharp constrast to the anisotropic distribution of light within 50 Mpc ||. 

There are two extreme explanations for the observed isotropy. On the one hand, it 
may happen that a bunch of sources that by pure chance are very close to us dominate 
the spectrum at the highest energies, and the particle orbits are bent [||]. This scenario 
requires large scale intervening magnetic fields with intensity 0(/j,G), to provide suffi- 
cient angular deflection. On the other hand, one can argue that there are many cosmic 
ray sources, even at the highest energies. The lack of plausible nearby sources in the 
arrival direction has encouraged the idea of positing undiscovered neutral hadrons, as 
well as mechanisms which are able to break the GZK barrier. Although sufficiently 
heavy particles would avoid the GZK cut off (the cut-off energy varies as the square 
of the mass of the first resonant state) ||, the existence of these particles now appears 
to be excluded by laboratory experiments [0]. The only standard model (SM) parti- 
cle that can reach our galaxy from high redshift sources without significance loss of 
energy is the neutrino. The expected event rate for early development of a neutrino 
shower, however, is down from that of an electromagnetic or hadronic interaction by 
six orders of magnitude. This problem can be raised by simply postulating that the 
total neutrino-proton cross section a vp increases at high center of mass energies, y/s > 
TeV. The hypothesis that all particles may have a strong interaction above collider 
energies is certainly not new ||. Recently, some scenarii with n large compact di- 
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mensions and precocious unification around the TeV-scale have rekindled this idea ||. 
Within this framework, SM fields are trapped into a 3+1 dimensional thin shell and 
only gravity propagates in the higher dimensional space. Therefore, the compactifi- 
cation radius r c of the extra dimensions can be large, corresponding to a small scale 
l/r c of new physics. Here, the weakness of gravitational interactions is a consequence 
of the large compactification radius, encoded in the relationship between the Newton 
constant G = M~ 2 and the fundamental scale of gravity M* ~ TeV, M 2 X ~ r™M™ +2 . 
From the 4-dimensional perspective, the higher dimensional graviton appears as an 
infinite tower of Kaluza-Klein (KK) excitations [|Kj. The weakness of the gravitational 
interaction can be thus compensated by the very large multiplicity of KK states. As a 
consequence, cross sections mediated by spin 2 particles increase rapidly with energy 
TTJ . The KK model may fail [12] in that the KK modes couple to neutral currents, and 



the scattered neutrino transfers only about 10% of its energy per interaction, thereby 
elongating the shower profile. Nevertheless, an s 2 growth of a up , supplemented by 
multiple scatters within the nucleus, yields enough energy transfer to save the model 

0- 

On a different track, if some flavor of neutrinos has masses m„. (9(10 _1 ) eV, the 1.9 K 
thermal neutrino background is a target for extremely high energy neutrinos to interact 
forming a Z-boson that subsequently decays producing a "local" flux of nucleons and 
photons [14 1 . The energy of the neutrino annihilating at the peak of the Z-pole is well 
above the GZK limit 

E„ = ^ = A f^lxlO-eV. (1) 
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The mean energies of the ~ 2 baryons and ~ 20 7-rays in each process can be estimated 
by distributing the resonant energy among the mean multiplicity of 30 secondaries. The 
proton energy is given by 

< E p > ~ -^L „ 1.3 (*L) x 10 20 eV, (2) 
60 m„, \m Uj J 



whereas the 7-ray energy is given by 

< E „ o.7 (—)x 10 20 eV. (3) 

7 120 77V, \m v J y ! 

The latter is a factor of 2 smaller to account for the photon origin in two body 7r° 
decay. The annihilation/Z-burst rate can be amplified if neutrinos are clustered rather 
than distributed uniformly throughout the universe |Tj|. In such a case the probability 
of neutrinos to annihilate within the GZK zone is on the order of 1%, with the exact 
value depending on unknown aspects of neutrino mixing and relic neutrino clustering 
(more on this below). 
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Adding to the puzzle, the AGASA experiment has already reported data strongly 
suggesting that the pairing of events on the celestial sky is occurring at higher than 
chance coincidence ||16|| . Specifically, four doublets and one triplet of showers with 
separation angle less than the angular resolution 2.5° are observed among the 36 events 
reported with mean energy above 4 x 10 19 eV. The chance probability of observing such 
a triplet in an isotropic distribution is about 1%. Another doublet is observed if we 
include events above 3.8 x 10 19 eV. The arrival directions in a combined data sample 
with three other surface experiments further suports non-chance association, especially 
in the direction of the SuperGalactic plane |l7j. If not a statistical fluctuation, the 
event clustering would have profound implications for models discussing the origin of 
ultrahigh energy cosmic rays. In this paper we elaborate on this issue. 



2 Estimates of cluster probabilities 

Let us begin by reviewing the current status of event clustering including the recently 
enlarged sample reported by the AGASA experiment fi~6fl . To proceed, we adopt the 
formalism introduced in Ref. |18|| . We start considering the solid angle Q on the 



celestial sphere covered by an experiment divided into N equal angular bins, each with 
solid angle u ~ tt8 2 . Then, by tossing n events randomly into 
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(4) 



bins, one is left with a random distribution. Now, we identify each event distribution 
by specifying the partition of n total events into a number m of empty bins, a number 
rrii of single hits, a number m 2 of double hits, etc., among the N angular bins that con- 
stitute the whole exposure. After a bit of algebra, it is easily seen that the probability 
to obtain a given event topology is 



3=0 



[rrii 



rrin 



(5) 



where 



~ = x '.v) TT- 



ie) 



Using Stirling's approximation for the factorials with the further assumption N ^> n ^> 
1, Eq. (||) can be re- written in a quasi- Poisson form 



P({mi},n,N) « V 



n 



\ m j) 3 -mTri(j-2)\ 



J=2 3 



(7) 
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where r = (N — m )/n « 1, and the prefactor V is given by 



For "sparse events" , where N 3> n, one expects the number of singlets m x to approxi- 
mate the number of events n. In such a case the prefactor is near unity. 

To estimate the celestial sky coverage Q for ground-based experiments, one rotates 
the fixed-coordinate solid angle about the earth's axis of rotation. For experiments 



with vertical acceptance from the zenith to 9 Z the relevant formula is |18] 

„_ J 2 sin 6 Z cos a, for a + 6 Z < 90° 

\l + sm(9 z -a), for a + 6 Z > 90° 

where a is the latitude of the experiment. 

The experimental data reported by Volcano Ranch, Haverah Park, Yakutsk and 
AGASA are used to determine the purely statistical probabilities for various cluster 
topologies. As recommended in Ref. W% only extensive air showers with zenith angle 
8 Z < 45° (that have good quality in the energy and arrival direction determination) 
are taken into account. The updated data list of possible cluster members is quoted 
in Table 1. Using Eqs. (|j) and ([5]) we calculate the inclusive probabilities for various 
cluster topologies as a function of the angular resolution and the accumulated number 
of events n. By inclusive probabilities we mean that the specified number of j-plets plus 
any other cluster, counts as all the j-plets + extra-clusters. The main experimental 
properties for the four experiments are summarized in Table 2. 

In Fig. 1 we show the inclusive probabilities for one and more doublets and two and 
more doublets for the sample of Haverah Park. The chance probability for clustering 
within 3° is larger than 50% and hence not statistically significant. In Fig. 2 we show 
the inclusive probabilities for 8 doublets and 2 triplets at different CR-sky coverages. 
The probability of chance association is only "small" for angular binning tighter than 
3°, and Q > 4. The chance probability for clustering within 4° and 5° remains always 
larger than 10%. Therefore, the observation of this topology within the approximate 
angular resolution of the combined data set, is not statistically significant. This result 



agrees with previous numerical simulations ||T7 . 

We now examine whether there is any evidence for clustering above the statistical 
expectation when considering the AGASA subsample. The latter has much better 
angular resolution. In Fig. 3 we show the chance probabilities of observing 5 doublets 
and one triplet given 58 events at AGASA. The probability is extremely sensitive to the 
angular binning. In this case, the chance probability within the experimental angular 
resolution is less than 10~ 3 . This result is in very good agreement with the one recently 
obtained using numerical simulations of the angular two point correlation function of 
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ultrahigh energy CRs: A 3 x probability of chance clustering with a bin size of 
2.5° and an energy cut-off at 4.8 x 10 19 eV |fL9| . Furthermore, when including data 
from Yakutsk's experiment above 2.4 x 10 19 eV, the combined probability of chance 
clustering is reported to be as small as 4 x 10~ 6 , strongly suggesting that CR sources 



are point-like on cosmological scales [19 



3 Correlation with high redshift objects 

Compact radio quasars (CRQSOs) are strong radio emitters, a fact that along with 
their variability, is indicative of strong beaming. The bulk of the observed non-thermal 
emission of these objects is thought to be produced in strong, relativistic jets of charged 
particles emitted by the active nucleus, which is likely formed by an accreting super- 
massive black hole. These powerful objects have been under suspicion as the primary 
source of ultrahigh energy CRs for some time now [20, 21, [22]. Therefore, we find it par- 



ticularly attractive to examine whether there exists a correlation between CR-clusters 
and CRQSOs. We shall use the 451 CRQSOs with flat spectrum and declination above 
— 10° degrees taken from the surveys of Ref. f23|. With the aim of finding the posi- 



tional coincidences and evaluating their significance, we adopt the procedure of Ref. 



2jfl . First, we look for real correlations between the two sets. In order to do so, we 
consider a circle around the centroid of each CR event, this circle has a radius equal 
to the reported 1 sigma positional error (see Table 2). If a CRQSO is within the circle 
of all members of the cluster, we say that there is a positional coincidence. We are not 
giving a higher significance to directional coincidences with small offsets than to coinci- 
dences that are not so close, just because the original errors of the CRs are of the order 
of degrees. As a first trial, we look for positional coincidences between the CRQSO 
sample and all the events listed in Table 1. We have found that there are no objects 
which correlate with the direction of triplets, so in the following a cluster denotes a 
pair of events. Just 4 clusters (out of 10) of the sample are positionally coincident 
with CRQSOs. Numerical simulations using large numbers of synthetic populations 
(thousands of them were made for each correlation study) sampled randomly and uni- 
formly in right ascension and declination, are then performed in order to determine 
the probability of pure chance spatial association. Strictly speaking, we generate syn- 
thetic populations of 451 CRQSOs and compare them with the actual positions of the 
CRs. We have taken into account that the artificial sets of CRQSOs are constrained 
(as are the actual ones) to the declination range 5 > —10°. In Fig. 4 we show the sky 
distribution of 451 CRQSOs in equatorial coordinates. The apparent anisotropy is due 
to obscuration from the galactic plane of our galaxy. Since only one ultrahigh energy 
CR cluster lies in this region, the real level of positional coincidences (even considering 
that an uniform distribution of CRQSOs is there) cannot significantly differ from what 
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is observed. The present Monte Carlo simulations preserve an isotropic distribution for 
background sources. The level of random positional coincidence after 5000 simulations 
(a larger number of simulations does not significantly modify the result) is shown in 
Fig. 5. Notice that the actual coincidences are less than 1 standard deviation away 
from the simulated mean value 3.5 ± 1.8. As a second trial, we repeat the whole anal- 
ysis but considering only CRs in the AGASA sample above 4 x 10 19 eV. In this case 
we found that the number of real matches is 3, whereas the expected number from 
pure chance estimated from the simulations is 0.45 ± 0.66, i.e., 3.8 standard deviations 
below (see Fig. 6). The Poisson probability of a random occurrence of any number of 
coincidences greater than or equal to the real positional coincidence is P = 9.7 x 10~ 3 . 



4 Outlook 

Let us end with a discussion (and some speculations) on the implications of our result. 
The clustering beyond statistics by itself imposes certain constraints on possible CR- 
sources. As can be read in Table 1, the event-pairing has members with rather different 
energies. The energy spread would have profound consequences for the propagation of 
charged ultrahigh energy CRs, even in the regular magnetic field of (9(nG). Strictly 
speaking, if the Larmor radius of a particle (r"L — 10 2 Mpc .E20/-B-9) is much larger 
than the coherence length of the magnetic field £ co h, the typical deflection angle from 
the direction of the source, located at a distance d, can be estimated assuming that 



the particle makes a random walk in the magnetic field |24 



where E20 is the energy of the particle in units of 10 20 eV, and _B_g is the magnetic 
field in units of 10~ 9 G. Now, it is straightforward to check that scatterings in large 
scale magnetic irregularities O (nG) |25| are enough to bend the orbits of trans-GZK 
protons more than 5° in a 50 Mpc trip. The time delay in the arrival of cosmic ray 
pairs can be used to impose additional constraints on the nature of the sources. In 
particular, the time lags shown in Table 1 are long enough to rule out an origin based 
on bursting sources. The average time delay in the pairs of the sample (with real 
positional coincidences) is 3 yr ±7 months, favoring a compact (but not transient) 
source. Typical source sizes for the production of high energy particles in AGNs are 
smaller that the radius of the outer gamma-spheres, i.e. < 10 15 cm [p6[| . All in all, if 
forthcoming data confirm that the clustering is not a statistical artifact, we are left with 
two main possibilities: (i) the particles are charged but then the sources should be very 
nearby (which is hard to accept already), or (ii) the CR constituents of the clusters are 
neutral particles! The first option is very problematic because if the sources are nearby 
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(say located within our own Galaxy) global isotropy should not be observed, contrary 
to the evidence, unless strong magentic fields could produce the required isotropization, 
in which case the clustering effect is broken. 

In order to explain the event-pairs (pointing towards distant sources) by neutrino 
showers with TeV-scale quantum gravity, neutrinos must be capable of generating 
vertical atmospheric cascade developments both at the lower and higher energy pair. 
The rapid energy behaviour of cross section in KK models makes this requirement 
rather difficult [0. The Z-burst, however, still appears as a viable model. As stated 
above the secondaries produced by Z decay are mainly photons. Because cascading by 
e + e~ pairs (produced by double pair production on the microwave background, or pair 
production on the radio background) is rapid, the photon energy is quickly reduced 
below the GZK-limit. Just to get an estimate, the mean interaction length at 3 x f 20 
eV is a few Mpc, and the energy attenuation length (which is very sensitive to the 
magnetic field because of synchrotron losses) assuming extragalactic magnetic fields O 
(nG) is ~ 20 Mpc. Both increase with decreasing energy. As mentioned above, some 
protons are able to survive a 60 Mpc trip with energies above 10 20 eV. Besides, the 
nucleon-channel would suffer magnetic field deflections reducing the correlation with 
the neutrino-emitting-source. For instance, _B_g ~ 2 would defocus a proton E p ~ 1.5 x 



10 20 eV (as the one detected by AGASA on 97/03/30) more than 5° f2|. Therefore 



since the photon flux is depleted through interactions with the radiation backgrounds, 
the relation 10:1 of photon/nucleon in the Z decay is significantly reduced, reducing 
also the correlation with CRQSOs (because the protons are bent in the magnetic field). 
In other words, the Z-burst model predicts a correlation between ultrahigh energy CRs 
and highly redshifted background sources ||22|| , which diminishes with increasing energy 
(say E > 8 x 10 19 at 1 a level) [|21| . We note that none of the members of clusters 
with real matches in our analysis satisfy the above energy cut-off. We do not attempt 
to make yet another estimate here and only note that on average, proton deflections 
of ~ 10° [28 are consistent with existing bounds on neutrino masses, relic neutrino 
density and magnetic field strengths. 

In closing we would like to note that the Z-burst model is not free of problems: 
Firstly, it is desirable to understand more completely the depletion of high energy 
photon flux, since the characteristcs of the air showers detected so far do not seem 



to resemble an electromagnetic 7-ray shower ||29|| . Secondly, there is the additional 
difficulty of getting neutrinos with energies ~ 10 22 eV. If these are secondaries from 
pion production, this implies that the primary protons which produce them must have 
energies of ~ 10 23 eV. While standard acceleration mechanisms require quite extreme 
parameters to achieve such extraordinary energy, we note that only a few dozen of such 
sources in the whole visible Universe would suffice. 

Mysterious: adj. Having an import not apparent to the senses nor obvious to the 
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intelligence; beyond ordinary understanding [BO . This appears to be the situation 
concerning the seeming ultrahigh energy cosmic ray clustering discussed in this paper. 
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Figure 1: Inclusive probabilities for various clusters in a 27 event sample as seen by 
the Haverah Park experiment. The solid line is the exact result, whereas the dashed 
line is the Poisson approximation. 
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Probability of 8 doublets + 2 triplets (103 events) 
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Figure 2: Inclusive probabilities for 8 doublets and 2 triplets in a 103 event sample for 
various celestial solid angles. 
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Figure 3: Inclusive probabilities for 5 doublets and 1 triplet in a 58 event sample at 
AGASA. Solid (exact), dashed (Poisson). 
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Figure 4: Sky distribution of 451 CRQSOs in equatorial coordinates. The sample is 
complete within the sensitivity of the radio surveys (see Ref. ||23|| ). The apparent 
anisotropy is due to obscuration from the galactic plane of our galaxy. 
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Figure 5: Simulated positional coincidence between all CR-clusters taken from Table 
1 and CRQSOs for 5000 runs. 

17 



i 1 1 1 r 



3000 , 



2500 , 



2000 , 



1500 , 




Coincidences between clusters and CRQSOs 



Figure 6: Simulated positional coincidence between AGASA-clusters with energies 
above 4 x 10 19 eV and CRQSOs for 5000 runs. 
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Cluster 


Experiment 


Date 


Log E 


R.A. 


Dec. 


Triplet #1 


Haverah Park 
AGASA 
AGASA 


810105 
931203 
951029 


19.99 
20.33 
19.71 


20.00 
18.91 
18.53 


20.00 
21.07 
20.03 


Triplet #2 


AGASA 
AGASA 
AGASA 


920801 
950126 
980404 


19.74 
19.89 
19.73 


172.30 
168.65 
168.44 


57.14 
57.58 
55.99 


Doublet #1 


AGASA 
AGASA 


910420 
940706 


19.64 
20.03 


284.90 
281.36 


47.79 
48.32 


Doublet #2 


AGASA 
AGASA 


860105 
951115 


19.74 
19.69 


69.03 
70.39 


30.15 
29.85 


Doublet #3 


Haverah Park 
AGASA 


860315 
960513 


19.71 
19.68 


267.00 
269.05 


77.00 
74.12 


Doublet #4 


Haverah Park 
Yakutsk 


720525 
911201 


19.65 
19.62 


239.00 
235.40 


79.00 
79.80 


Doublet #5 


Volcano Ranch 
Haverah Park 


610319 
850313 


19.73 
19.62 


154.10 
157.00 


66.70 
65.00 


Doublet #6 


Haverah Park 
Yakutsk 


661008 
750317 


19.67 
19.67 


164.00 
163.70 


50.00 
52.90 


Doublet #7 


Haverah Park 
AGASA 


740228 
980330 


19.86 
19.84 


264.00 
259.16 


58.00 
56.32 


Doublet #8 


Haverah Park 
Haverah Park 


760206 
850313 


19.62 
19.62 


165.00 
157.00 


64.00 
65.00 


Doublet #9 


AGASA 
AGASA 


960111 
970410 


20.16 
19.58 


241.5 
239.50 


23.00 
23.70 


Doublet #10 


AGASA 
AGASA 


961224 
000526 


19.70 
19.70 


213.75 
212.00 


37.7 
37.1 



Table 1: Updated list of triplets and doublets within space angles of 3°, 4° and 5°. 



Experiment 


Begin 


End/status 


Latitude 


Longitude 


Q (sr) 


^min 


n 


AGASA 


1990 


in operation 


35° 47' N 


138° 30' E 


4.8 


1.8° 


58 


Haverah Park 


1968 


1987 


53° 58' N 


1° 38' W 


3.7 


3.0° 


27 


Yakutsk 


1972 


in operation 


61° 42' N 


129° 24' E 


3.1 


3.0° 


12 


Volcano Ranch 


1959 
1972 


1963 
1974 


35° 09' N 


106° 47' W 


5.0 


3.0° 


6 



Table 2: Event rates, celestial solid angles Q (sr), and angular resolutions # m in- We 
remark that the celestial solid angles have been reduced by a factor of 0.7 (see j76|/ for 
details). 
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